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ABSTRACT 


This  thesis  is  concerned  with  solutions  of  a  moving  boundary 

problem. 


A  few  proofs  of  existence  and  unicity  of  solution  to  the  problem 
are  outlined. 

A  formal  solution  due  to  Goursat  which  uses  Holmgren's  lemma 
is  presented. 

Several  analytic  and  non=analytic  techniques  which  may  ba  used 
to  obtain  numerical  solutions  are  presented. 
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CHAPTER  I 


INTRODUCTION 

In  this  chapter  we  state  the  problem  in  general  terms  and 
try  to  classify  the  physical  problems  into  well  defined  types.  We 
also  note  a  few  interesting  aspects  of  the  problem.  We  give  a  brief 
survey  of  the  remaining  chapters  and  explain  some  of  the  notations 
used  later  in  the  thesis, 

1.1  Introduction 

There  are  a  number  of  problems  in  heat  conduction  in  which 
one  material  is  transformed  into  another  phase  with  the  generation  or 
absorption  of  heat.  Examples  of  such  phase  change  processes  are  melt¬ 
ing,  evaporation,  recrystallization  and  dissolution.  In  these  processes 
one  often  encounters  a  boundary  value  problem  of  the  heat  equation  in 
a  domain  with  a  boundary  which  is  not  fully  known  and  whose  determina¬ 
tion  is  an  essential  part  of  the  solution.  Such  boundaries  are  referred 
to  in  the  literature  as  free  or  moving  boundaries.  We  will  refer  to 
the  boundary  as  moving  to  distinguish  it  from  similar  boundaries 
encountered  in  hydrodynamics  where  one  encounters  a  free  boundary 
between  dead  and  moving  liquid  surrounding  a  moving  body  in  water. 

In  phase  change  problems  the  boundary  between  the  material  in  the  two 
phases  moves  as  some  of  the  material  changes  from  one  phase  to  the 


other . 
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With  a  mewling  boundary,  the  problem  becomes  nonlinear  but 
that  is  not  the  only  complication.  The  problem  is  further  complicated 
by  the  fact  that  the  thermal  properties  of  the  materials  in  the  two 
phases  are  not  usually  the  same. 

The  governing  equation  for  the  temperature  distribution, 
u(x,t),  x  denoting  position  and  t  time,  is  parabolic.  Formal 
solutions  of  the  problem  depend  almost  entirely  on  a  lemma  due  to 
Holmgren  (see  Holmgren  (1908)^,  Goursat  (1915)>  Kolodner  (1956)). 
Holmgren’s  lemma  has  apparently  been  overlooked  in  much  of  the  recent 
literature.  We  note  that  formal  solutions  cannot,  in  general,  be 
expressed  in  terms  of  computed  functions.  Thus  we  resort  to  numerical 
computation  of  the  solution.  In  developing  a  numerical  solution  we 
need  an  appreciation  of  the  behavior  of  the  desired  function  on  and 
near  the  moving  boundary, 

1.2  The  mathematical  problem 

We  can  state  in  general  terms  the  mathematical  nature  of  the 
problem  considered.  Referring  to  Figure  1.2.1,  we  seek  a  function 
u(x, t)  in  regions  I  and  II  ,  Region  I  is  bounded  by  the  lines 
t  =  t^  and  t  =  tg  and  by  a  fixed  plane,  say  x  =  0,  and  by  the 
curve  T  whose  equation,  in  general  unknown,  is  x  =  X(t).  Region  II 
is  bounded  by  the  same  lines,  by  f  and  by  a  second  fixed  plane  which 


The  date  in  parentheses  following  the  name  of  an  author  indicates 
the  year  in  which  the  work  to  which  we  refer  was  published.  The 
work  is  found  in  the  bibliography  by  first  looking  for  the  author 
and  then  looking  for  a  work  published  in  the  year  given. 
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may  be  at  a  finite  or  infinite  distance.  The  governing  equations  for 
u  in  the  two  regions  are 


(1.2.1) 


c,u 
1  XX 


U. 


c„u  =  u 
2  xx  t 


The  numbers  and  are  physical  constants  which  depend  on  the 

material.  The  function  u  is  given  on  the  boundary  F  .  The  derivative 

u  may  be  required  to  be  continuous  or  to  possess  a  discontinuity  of  a 
X 

specified  amount  (e.g.  if  latent  heat  is  given  up)  on  crossing  F  . 

The  remaining  initial  and  boundary  conditions  are  familiar  and  irrele¬ 
vant  to  the  present  exposition^  they  will  be  taken  up  later. 


The  essential  difficulty  in  the  problem  is  the  determination 
of  the  unknown  moving  boundary,  X(t)  .  The  specification  of  u  on  F 
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is  a  non-linear  boundary  condition  and  the  problem  is  non-linear  even 
though  the  differential  equations  are  linear.  Let  us  consider  the 
matter  of  non-linearity.  Suppose  we  have  the  following  condition 

(1.2.2)  -u  (X(t),t)  =  H(t)  . 

It  is  clear  that  the  problem  is  non-linear  since  if  we  have  two  different 
functions  H^(t)  and  Hg(t),  we  obtain  two  different  boundaries  X^(t) 
and  X  (t)  .  The  two  solutions  would  not  apply  in  the  same  region  of 
the  (x,t) -plane  and  the  solution  for  the  sum  H^(t)  +  Hg(t)  cannot 
be  obtained  by  superposition, 

1.3  Classification  of  physical  problems 

A  classification  of  various  possible  physical  problems  will 
be  helpful  in  later  considerations.  The  problems  may  be  classified  by 
the  initial  existence  or  non-existence  of  regions  I  and/or  II.  We 
consider  six  types  and  give  examples.  The  classification  is  probably 
not  complete  but  it  gives  physical  content  to  the  mathematical  state¬ 
ment  of  the  problem. 

1.3»1  Only  region  I  exists  and  the  equation  of  the  moving  boundary 
may  or  may  not  be  given.  The  building  of  a  clay  embankment  provides 
two  problems.  If  the  rate  of  building  is  known,  thlen  the  problem  is  to 
determine  the  pore  pressure.  Alternatively,  the  problem  may  be  to 
determine  the  rate  of  building  such  that  the  pore  pressure  does  not 
exceed  a  specified  value. 
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1.3*2  Both  Regions  I  and  II  exist  initially,  the  inital  states  being 
defined.  Such  a  problem  is  exemplified  by  the  introduction  of  a  large 
mass  of  supercooled  ice  into  a  large  mass  of  water.  Light foot  (1929) 
gives  solutions  by  appealing  to  physical  reasoning. 

1.3*3  Region  I  does  not  exist  initially  but  is  continuously  created 
at  the  expense  of  Region  II,  A  classical  example  here  is  the  depth  of 
ice  in  a  permafrost  region  (Carslaw  and  Jaeger  (1959').  Other  examples 
are  the  formation  of  ice  on  cylindrical  cables  and  the  freezing  of 
ingots  (Miles  (1950)). 

1*3*4  Region  I  does  not  exist  initially  and  is  destroyed  as  it  is 
created.  Landau  (1950)  discusses  problems  of  this  type.  An  example 
is  the  melting  of  a  mass  of  metal  where  a  heat  source  is  applied  to 
one  face  and  the  molten  metal  is  continuously  blown  away. 

1.3*5  Region  I  is  continuously  created  at  the  expense  of  Region  II 
which  is  maintained  in  a  uniform  and  homogeneous  state.  Problems  of 
this  type  are  approximations  to  those  of  type  1*3*3  above.  In  the 
literature,  these  problems  are  referred  to  as  Stefan  problems. 

1,3*6  Neither  Region  I  nor  II  exist  initially  and  the  initial  position 
of  the  moving  boundary  is  not  known.  This  type  of  problem  is  found  in 
consolidation  theory.  Since  the  pore  water  moves  with  finite  velocity, 
the  initial  state  is  neutral;  consolidation  or  swelling  will  exist  at 
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t+0  but  not  at  t  =  0  .  Likewise,  the  boundary  at  t  =  0  is  not 
defined  but  its  limiting  position  as  t~»0  can  be  determined. 

1.4  Survey  of  thesis 

The  problem  we  consider  has  a  few  interesting  aspects.  Before 
we  can  examine  any  solutions,  we  examine  the  conditions  under  which  a 
solution  exists  and  is  unique.  After  doing  so  in  the  first  part  of 
Chapter  2,  we  obtain  a  formal  solution  which  requires  the  application 
of  Holmgren's  lemma  for  the  evaluation  of  the  discontinuity  at  the 
moving  boundary.  In  chapter  3  we  consider  several  techniques,  both 
analytic  and  non-analytic  and  having  varying  ranges  of  application. 

Some  of  these  techniques  are  those  which  belong  to  the  hand  calculation 
stage  in  the  solution  of  the  problem.  A  few  techniques  may  be  used  as 
machine  techniques. 

1.5  Notation 

Before  proceeding  we  explain  the  notation  we  will  use.  Sub¬ 
scripts  t  or  x  on  any  function  indicate  differentiation  with  respect 
to  that  variable.  In  general  X(t)  will  be  used  for  the  moving  boundary. 

We  will  use  the  following  function 

(exp  (-{x-r )”)/( 4(t-s ) ) )  /(s/t>s) 

frequently.  For  convenience  we  will  denote  this  by  E(x,tjr,s). 


D.  E.  and  A.  E.  are  used  to  mean  differential  equation* 


}  lit  j  5s, 
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and  difference  equation  respectively. 

We  use  erf  (x),  the  usual  abbreviation  for 


erf  (x)  = 


,x 


exp 


dt 


) 


the  error  function.  We  use  erfc  (x)  for  1  ~  erf  (x),  the  complement 
of  the  error  function.  We  use  ierfe  (x)  to  represent 


ierfe  (x)  =  /  erfc  (t)  dt 

J  x 
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CHAPTER  2 


EXISTENCE  AND  UNICITY  OF  SOLUTION  AND  A  FORMAL  SOLUTION 

In  this  chapter  we  examine  the  conditions  under  which  the 
problem  has  a  solution  which  is  unique  and  indicate  the  methods  of 
proof  given  by  several  authors.  We  derive  a  formal  solution  to  the 
problem  giving  particular  attention  to  the  behavior  of  the  solution 
on  the  unknown  moving  boundary. 

2.1  Existence  and  unicity  of  solution 

Several  authors  have  given  proof  of  existence  and  unicity 
for  various  forms  of  the  moving  boundary  problem.  In  each  case  we 
state  the  problem  and  give  an  indication  as  to  the  method  of  proof, 

2.1.1  Evans'  proof 

Evans  (1951)  gives  a  proof  for 


(2. 1.1.2) 


u(x(t),t)  =  0 


(2. 1.1. 3) 


u 


X 


(0,t)  = 


(2. 1.1. 4) 


(2.1. 1.5) 


X(0)  =  0  . 


i  ■  :v  c  . 
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Evans  requires  that  u  have  bounded  continuous  first  derivatives  for 
0  <  x  <  X{t ) . 


In  the  proof 


u  )  dxdt  =  0 
t 


is  evaluated  by  application  of  the  boundary  conditions. 
u(x,t)  are  determined  iteratively  by  setting 


(t) 


u 


n~l 


(x,t)  dx 


X(t)  and 


which  satisfies  (2. 1.1.4)  and  in  the  limit 

lim  u  ^(x,t)  =  u(x,t)  . 
n  -*  oo 


Unicity  is  established  by  assuming  the  existence  of  two 
solutions  for  the  boundary,  X(t),  Y(t),  X(t)  >  Y(t)  .  Thus 


X(t)  =  t 


Y(t)  =  t 


u(x,t)  dx 


v(x,t)  dx 


We  form  the  difference 


rY(t)  rX(t) 

X(t )  -  Y(t)  =  /  (v(x,t)  -  u(x,t))  dx  =  /  u(x,t)  dx 

-'r,  J  Y(t) 


We  let 
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r  %  **  , 
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w(x,t)  =  v(x,t)  -  u(x,t) 
which  satisfies  w  (x,t)  =  w  (x,  t) 

XX  L 

with  boundary  conditions 

wx(0,t)  =  0  and  w(Y(t),t)  =  -u(Y(t),t) 

and  -u(Y(t),t)  <  0  since  u(x,  t)  >  0  . 

Since  wx(x, t)  =  0,  we  can  reflect  this  solution  about  x  =  0  so  that 
w(x,t)  still  satisfies  the  D.E.  but  with  boundary  conditions 

w(-Y{t),t)  =  -u(Y(t),t)  <  0 

and 

w(Y(t),t)  =  -u(Y(t),t)  <  0 

By  the  maximum  principle  for  parabolic  equations. 

If  u  is  a  solution  of 

A(x,t)u  -  B(x,t)u  -  C(x,t)  =  0 
where  the  coefficients  are  bounded  and  continuous, 

A,  B  >  0,  then  G  <  0  (C  >  0)  implies  u  attains 
its  positive  maximum  (negative  minimum)  on  the 
boundary  of  the  domain, 

w(x,t)  <  0  since  A  =  1,  B  =  1,  €  =  0  . 

Thus 

rY(t)  rX(t) 

x(t)  -  Y(t)  =  /  w (x, t )  dx  -  /  u(x,t)  dx  <  0 

Jo  '  ^Y(t) 
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which  contradicts  the  assumption  X(t)  >  Y(t)  ,  Hence  X(t)  =  Y(t)  . 


2.1.2  Kolodner' s  proof 

Kolodner  (1956)  considers  the  problem 

(2. 1.2.1)  uxx(x,t)  =  ut(x,t) 

(2. 1.2.2)  u(X(t),t)  =  f(t)  ,  |f(t)  |  <  A 

(2. 1.2. 3)  ux(x(t),t)  =  g(t),  |g(t)|  <  At*1/2 

Thus  f  and  g  are  functionals  of  X(t),  f  being  continuously 
differentiable  and  g  Lipschitz  continous. 

Kolodner  establishes  unicity  before  going  on  to  show  existence. 
In  so  doing  he  shows  that  if  (2. 1.2. 2)  or  (2. 1.2. 3)  has  at  most  one 
solution  then  the  moving  boundary  problem  has  at  most  one  solution.  He 
assumes  two  solutions,  u^,  u^  and  sets  u  =  u^  -  u^,  u  also  being 
a  solution.  By  employing  the  identity 


// 


f  dxdt  = 
x 


Domain 


/ 

Boundary 


(f  dx  +  2ff  dt) 

X 


he  shows  u  =  0  and  thus  . 


In  proving  existence  he  shows  that  the  problem  can  be  reduced 


to  a  functional  equation  for  the  moving  boundary  and  shows  that  the 


12 


existence  of  a  solution  X(t)  insures  the  existence  of  a  solution  to 
the  moving  boundary  problem  for  X  <  0  . 

Looking  at  the  problems  considered  by  Evans  and  Kolodner  it 
would  appear  that  Kolodner  treats  a  more  general  problem. 

2.1.3  Douglas*  proof  of  unicity 

Douglas  (1957)  establishes  unicity  for 


(2. 1.3.1) 

u^fot)  »  f(w)  ut(x,t)  , 

f(u)  >  0  , 

(2. 1.3. 2) 

Ux(0,t)  m  -h(t) 

(2. 1.3.3) 

u(X(t),t)  =  0 

(2. 1.3. 4) 

0 

11 

0 

(2. 1.3. 5) 

Xfc(t)  =  b  -  cux(X(t),t)  , 

b,  c  >  0  , 

He  uses  the  maximum  principle  for  parabolic  «quations.  In  our  problem 
u  attains  its  maximum  (minimum)  on  the  boundary. 

2.1.4  Mir anker's  proof 

Miranker  (1958)  shows  that  the  following  problem  has  a  solution 
uxx(x,t)  =  ut(x'fc) 
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ux(0,t)  «  -g(t)  ,  g(t)  >  0 

u(X(t),t)  =  0 

u{x,0)  =  f(x)  ,  f(x)  >  0 

ux(x(t),t)  =  -Xt(t)  , 

where  g,  f  e  C  ,  g  being  differentiable  and  f  continuously  different 
tiable  for  0  <  x  <  A  ,  f(x)  =  f(-x)  ,  f(x)  =  0  for  x  >  A  } 
f  (A)  <  0  ,  f(0)  =  -g(0)  . 

A 

An  example  of  an  f  satisfying  the  above  is 

f(x)  =  g(0)x2/A  -  g(0) |x |  +  g(0)A  -  A2  ,  Jx |  <  A 
f  (x)  =  0  ,  |x  |  >  A  . 

Miranker  uses  Kolodner's  method  to  obtain  a  functional  equation  and  then 
uses  the  principle  of  contraction  mapping  to  solve  the  problem. 

2.1.5  Kyner 1 s  work 

Kyner  treats  the  problem 

u„(x,t)  =  f(x,t,u)  ut(x,t) 

ux(0,t)  =  -g{t)  , 

ux(x(t)# t )  =  -  xt(t) 


g(t)  >  0  ,  g  e  C 
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X(o)  =  A 

u(x,0)  =  h(x)  ,  h  e  C  ,  0  <  f(x)  <  b(A-x)  ,  b  >  0  * 


He  establishes  unicity  in  a  manner  similar  to  that  of  Douglas  (1957) • 
To  prove  existence,  he  uses  the  relation 


g(t)  =  c(t) 


F(x, t ,u(x, t,X(t ) ) )  dx 
X(s ) 

H(x,s,u(x,s,h))  dxds 


ru 

where  F(x,t,u)  =  /  f(x,t,s)  ds 

J  o 


and  H(x,t,u)  =  -  /  ffc(x,t,s)  ds  , 

Jo 


to  define  a  mapping  of  a  set  of  boundary  curves  into  itself.  Then  there 
is  at  least  one  boundary  curve  which  is  invariant  under  the  mapping  and 
which  gives  the  solution. 


2.1.6  Summary 

The  conditions  under  which  Kolodner  gives  a  proof  of  existence 
and  unicity  appear  to  be  the  most  general  and  most  widely  applicable. 

The  other  authors  treat  various  forms  of  a  Stefan  problem  which  is  a 
particular  type  of  moving  boundary  problem. 


2.2  A  formal  solution 

Having  considered  existence  and  unicity  of  solution  of  our 
problem,  we  now  consider  a  formal  solution.  The  approach  is  that  found 


_  j 
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in  Goursat  (1915)*  By  exhibiting  a  solution,  existence  of  a  solution 
is  established.  Unicity  follows  by  assuming  the  existence  of  two  such 
solutions  and  obtaining  a  contradiction. 

In  order  to  illustrate  some  of  the  difficulties  encountered 
in  solving  the  moving  boundary  problem  we  outline  a  formal  solution  of 
the  problem  paying  particular  attention  to  the  behavior  of  the  solution 
on  and  near  the  boundary  and  to  the  type  of  discontinuity  inherent  in 
the  solution  on  crossing  the  boundary. 

Before  proceeding  we  require  definitions  of  regular,  analytic 
and  holomorphic  functions  of  a  complex  variable.  A  function  is  regular 
in  a  region  if  it  is  continuous  and  has  continuous  first  partial  deri¬ 
vatives  in  that  region.  A  function  is  analytic  in  a  region  if  it  is  a 
regular  function  which  satisfies  the  Cauchy -Riemann  equations  in  the 
region;  i.e.,  if  f(x)  =  P(x,y)  +  i  Q(x,y)  then 

U  =  &  and  £  _  . 

dx  ay  oy  ox 

A  function  f (z),  is  holomorphic  in  a  convex  region  A  if,  in  addition 
to  being  analytic  in  A,  f(z)  is  such  that  to  each  z  in  A  there  corre¬ 
spond  unique  values  of  f(z)  and  f'(z)  . 

2.2.1  A  fundamental  solution 

We  consider  a  formal  solution  of 

uX3^X,t)  =  ut(x^t)  • 


(2.2.1. 1) 


:  t* 
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Suppose  we  make  the  transformation 


(2. 2. 1.2) 


Performing  the  calculations  we  find  that  (2. 2.1*l)  is  replaced  by  an 
equation  of  the  same  form 


(2. 2. 1.3) 


u '  =  u'  . 

XX  t 


Thus,  if  u(x,t)  is  a  solution  of  (2.2.1. l)  then  so  also  is 


Let  us  now  consider  an  integral  depending  on  two  arbitrary 
parameters  r  and  s  , 


E(x,t;  r,s) 


which  Goursat  (1915)  remarks  plays  a  role  relative  to  (2.2. l.l),  a  par¬ 
abolic  equation,  analogous  to  that  of  log  r  in  the  theory  of  harmonic 
functions . 

We  note  that  E(x,t;  r,s)  is  real  valued  for  t  >  s  .  Defin¬ 
ing  U(x,t;  r,s)  as  follows 


U(x,t;  r,s)  =  E(x,t;  r,s)  t  >  s 


0 


t  <  s 


. 

:  : 


' 
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we  see  that  it  is  a  solution  of  (2. 2.1.1)  which  is  continuous  and  has 
continuous  first  partial  derivatives  for  all  real  values  of  x  and  t 
except  x  =  r,  t  =  s  .  If  we  differentiate  U  with  respect  to 
x,  t,  r,  s  any  number  of  times  we  obtain  which  has  the  same 

properties  as  U  .  Each  of  the  is  a  solution  of  (2. 2.1.1) 

which  is  continuous  and  has  continuous  first  partial  derivatives  in 
the  entire  plane  except  at  x  =  r,  t  =  s  .  We  call  U  the  funda¬ 
mental  solution. 

From  the  fundamental  solution  we  can  derive  new  solutions 
by  integration.  Let  P(r,s)  ,  Q(r,s)  be  any  two  continuous 
functions  of  the  coordinates  (r,s)  of  a  point  M  along  an  arc  of 
a  curve  G  ,  r,  s  <  »  .  The  line  integral 

v(x,t)  =  J'  U(x,t;  r,s)  (P(r,s)  dr  +  Q(r,s)  ds) 
c 

is  a  solution  of  (2.2. l.l)  which  is  continuous  and  has  continuous  first 
partial  derivatives  in  every  domain  D  having  no  point  in  common  with 
the  curve  C  .  From  the  definition  of  U  we  take  for  the  calculation 
of  v(x,t)  only  that  portion  of  C  lying  below  the  line  t  =  t^  . 

If  t  =  t^  passes  through  the  minimum  ordinate  of  C  ,  then  v(x, t)  =  0 
Likewise  for  all  points  of  t  =  t^  not  passing  through  C  .  For  some 
line  t  =  t 1  where  C  is  below  t  =  t5  ,  v(x, t)  is  a  holomorphic 


function  of  the  two  variables  x  and  t 


■  vW 

M  . 

* 

■■  ;k  08  ws  .  :;:j:  -  , :  , 

- 
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2.2.2  Fundamental  solution  on  the  boundary 

The  most  important  aspect  of  our  problem  is  the  moving 
boundary  and  the  behavior  of  any  solution  on  or  near  the  boundary.  In 
this  section  we  consider  the  behavior  of  integrals  of  the  type 

(2.2.2. 1)  F(x,t)  =  [  f(t)  exp(-(x-X(t))2/4(t-r))/(t-r  j^2  dt 

J  a 

=  r  f(t)  E(x,t;  X(t),r)  dt 
*'  a 

where  f(t)  and  X(t)  are  continuous  on  the  closed  interval  [a,b] 
and  X(t)  represents  the  arc  C  .  F(x,t)  is  a  solution  which  is 
continuous  along  with  its  partial  derivatives  of  all  orders  in  every 
region  of  the  plane  having  no  point  in  common  with  X(t)  .  F  is 
continuous  in  x  and  t  on  X(t)  since  it  is  uniformly  convergent 
in  the  neighborhood  of  any  point  of  X(t)  , 


Suppose  we  evaluate  F  along  DE,  an  arbitrarily  small 
arc  of  X(t)  .  Then  |F|  <  e  independent  of  the  position  of  a 
neighboring  point  (x,y)  .  If  |f(r)  |  <  M  ,  then  along  BE, 


|F(x, t)  |  <  M 


<  e’ 


for  |t-a|  <  5  . 

crossing  X(t)  . 


We  observe  later  that  F  is  discontinuous  on 

x 


The  lines  t  =  t^,  t  =  t^  passing  through  the  lowest  and 
the  highest  point  of  X(t)  divide  the  plane  into  several  regions. 


' 


! 
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Below  and  on  t  =  ,  F(x,t)  =  0  ;  above  t  =  ,  F(x,t)  is  a 

holomorphic  analytic  function  of  x  and  t  .  We  are  interested  in  the 
nature  of  F  between  t  =  t^  and  t  =  t^  and  to  the  right  and  left 
of  X(t)~. 


Suppose  (x,t)  moves  on  x  =  a  .  F(a,t)  and  its  deri¬ 
vatives  with  respect  to  t  are  all  continuous  functions  but  it  is  not 
in  general  an  analytic  function  of  t  for  f(r)  any  continuous  function. 
The  value  of  F  along  x  =  a  depends  only  on  values  of  f(r)  along 
x(t )  for  t  <  t^  .  If  we  replace  f  by  another  continuous  function 
f^  such  that 

f(r)  =  f^(r)  along  X(t),  t  <  , 

4  f2_Cr)  along  X(t),  t  >  t^  , 

then  F(x,t)  will  coincide  with  F(a,t)  for  t  <  t  but  will  differ 

1 


(  c  ■■ 

'  "  ’•*  t'-  v  •  ;t:  .  /  :  ■  ;  •:  r 
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from  F(a,t)  for  t  >  t^  .  Since  is  arbitrary,  F(a,t)  cannot 

be  a  holomorphic  function  of  t  along  x  =  a  nor  even  on  any  portion 
of  x  =  a  if  f(r)  is  any  continuous  function.  However,  the  derivatives 
of  F(a,t)  satisfy  certain  inequalities  analogous  to  those  which  charac¬ 
terize  analytic  functions.  If  we  let 


(2. 2. 2. 2) 


G{x,t ) 


t  * 


X(r),r)0 


t-r  } 


dr 


the  F  =  -G  /2 

x  x' 


and  G  is  also  a  solution  of  (2. 2.1.1)  which  is  continuous  and  has 
continuous  first  partial  derivatives  for  a  <  t  <  b  except  on  X(t)  , 
the  path  of  integration.  We  extend  the  definition  of  F  to  the  entire 
plane  by  taking 

f(r)  =  0  for  r  >  b  . 

Thus  G  is  analytic  in  x  for  fixed  t  .  From  the  derivation  of 
F(x,t)  ,  the  derivatives  of  F(x,t)  are  majorized  by 


dxn 


<  Mnl/dn 


at  every  point  of  x  =  a  where  jF(x,t)|  <  M  for  | x-a  j  <  d  . 

G(x, t)  has  a  discontinuity  along  X(t)  which  Holmgren  has 
Assuming  that  Xt(t)  is  continuous  on  [a,b]  we  consider 


studied . 
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(2. 2.2. 3) 


H(x,t)  =  f  E(x,t;  x(r),r) 


ftrlLll) 

v  t-r  ' 


dr  . 


Since  it  does  not  have  a  simple  geometric  interpretation  we  study  its 
behavior  directly.  If  (x,t)  =  (X,T)  is  on  X(t)  ,  H(x,t)  has  a 
meaning 

H(x,l)  =  f  E(x,r;  X(r),r)  (x~^r^)  dr  . 

J  a 


However,  H(x, t)  approaches  different  limits  when  (x, t)  approaches 
(X,T)  depending  on  whether  (x,t)  is  to  the  right  or  left  of  X{t)  . 
To  show  this  we  consider 


(2. 2. 2.U) 


-2X’(r)  E(x,t;  X(r ),r )  dr 


which  is  continuous  on  X(t)  .  We  form 


H(x, t)  +  H^(x,t)  =  4 


exp  (-U  )  du 


where  u  = 

position  of  (x,t)  .  If 
When  r  varies  from  a 


and  c  and  d  are  determined  according  to  the 


x  >  X(t)  then  (x,t) 

to  t,  u  varies  from  x. 

1 


is  to  the  right  of  X(t) 


If  x  <  X(t)  ,  then  d  =  -c»  ,  Thus  for  (x, t)  not  on  X(t)  , 


r± 00  2 

[x, t)  +  H^(x, t)  =4  /  exp  (“U  )  du 

J  x. 


-  ' 


.3  ’  3  ’i  auriT 


-  22 


which  when  we  let  (x, t)  approach  (X,T)  becomes 


p±  00  2 

(2.2.2. 5)  lim  H(x,t)  +  lim  H^x^t)  =  4  /  exp  (-u  ) 

J  x 


du 


the  +(-)  taken  for  (x,t)  to  the  right  (left)  of  X(t)  .  If 

and  0 


x  =  X,  t  =  T,  the  limits  for  u  are  x„  * 


Thus 


(T-a 


o  o  o 

(2.2.2 .6)  H(X,T)  +  HX(X,T)  =4  /  exp  (-u  )  du  . 

J  V 

2 

If  we  subtract  (2. 2. 2. 5)  and  (2. 2. 2. 6)  term  by  term  and  note  that 
H^(x, t)  is  continuous  at  (X,T),  we  get 

p±  00  2 

(2. 2. 2. 7)  lim  F(x,t)  =  F(X,T)  +  4  /  exp  (-u  )  du 

■'  o 

=  F(X,T)  +  2(«)1//2  , 

the  +(-)  taken  when  (x, t)  is  to  the  right  (left)  of  x(t). 

Having  established  (2. 2. 2. 7)  we  consider  another  form  of 
G(x,t)  written  as  follows: 

(2. 2. 2. 8)  G(x  t)  =  ^  (f(r)“f(t))(x-X(r))  E(x,t]  X(r),r) 

J a  (*>*) 


fft)  rfc  (x-X(r))  E(Xjt;  X_(r ) ,r ) 

'  ;  J  a  (t-r) 


dr 


+  I2  • 


dr 


a  ■  2) 

■  V  .  .  ...  '  ,  .  ,  . 
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If  f(t)  satisfies  a  Lipschitz  condition  in  [a,b]  ,  then 


|ixl  <  M  / 

J  a 


dr 


/(t-r)1/2 


Consequently  is  uniformly  convergent  in  any  neighborhood  of  any 

point  (X,T)  on  X(t)  and  represents  a  continuous  function  in  the 
neighborhood  of  X(t)  .  The  discontinuity  in  G(x, t)  as  (x,t) 
approaches  (X,T)  is  inferred  from  the  discontinuity  of  Ig  . 

Thus  we  have  the  general  formula 

lim  G(x,t)  =  +  2jk  f  (l)  +  f  Xir) iT)  dr 

J  a  (T“r) 

where  +  or  -  is  chosen  as  above . 

If  X(t)  is  a  segment  of  x  =  a  ,  then  I0  =  0  and 


u(x  t)  -  —  f'  f(r)(x-a)  E(x,t;  a,r)  dr 

1  ’  J  af.  J a  (t-t) 


is  +  f(t)  when  (x,t)  approaches  (a,T)  on  x  =  a  . 

Holmgren’s  results  have  been  extended  by  Kolodner  (1956)  who 
observed  that  certain  of  the  conditions  on  X  and  f  are  only  required 
for  (0,t]  rather  than  [0,t].  Writing 

G(X,t)  m  r  +  f 

J  O  J  6 


. ,  {*)?  H> 

.  .  -3dfi  sfi  x-.QSQdo  ai  f,  ;  . 

..  I  £.'  .  - 
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1(0,6)  +  1(6, t) 


and  H(x,t) 


t 


1^0,6)  +  1^6, t)  , 


Kolodner  observed  that  both  1(0,6)  and  I^(0,6)  are  solutions  analytic 
in  x  and  t  for  t  >  6  >  0  while  in  [6,t]  Holmgren’s  conditions  are 
satisfied.  Holmgren  requires  that  f  be  a  function  such  that  |f|  <  A  and 
f(t)  is  Lipschitz  continuous  and  X  is  a  function  such  that  Xfc(t)  is 
continuous  in  the  closed  interval  [0,t]  whereas  Kolodner  requires  f  and 
X  which  satisfy  the  above  conditions  in  the  open  interval  (0,t]  . 

Thus  we  have  obtained  a  formal  solution  for  the  moving  boundary 

problem 


u(x,t) 


dr 


a 


which  is  an  integral  equation.  This  result  focuses  attention  on  the 
behavior  of  the  desired  function  near  a  moving  boundary  and  on  the  dis¬ 
continuity  inherent  in  the  solution  of  a  parabolic  equation  on  crossing 
a  boundary.  The  theoretical  importance  of  Holmgren’s  result  and  the 
extension  due  to  Kolodner,  and  of  F  and  G  in  the  analysis  of  boundary 
value  problems  is  undoubted.  The  practical  importance  of  F  and  G 
can  be  established  only  by  showing  that  they  can  be  used  to  determine 
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computable  solutions  in  applications  of  the  theory.  In  the  next  chapter 
we  consider  both  analytic  and  non-analytic  techniques  used  to  obtain 
solutions  for  a  number  of  problems. 


j.;.  J  rorij  ;r.?.  ■  :  v.,  ■  >sjrl  -  t-  x  ;>  s  ;h:.5  ’mfm ; 
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CHAPTER  3 

SURVEY  OF  PROBLEMS  AND  METHODS  OF  SOLUTION 

In  this  chapter  we  consider  both  analytic  and  non-analytic 
techniques  for  the  solution  of  problems.  The”  analytic  techniques  include: 

(1)  X(t)  a  t^7^  ,  (2)  Laplace  transforms  (3)  power  series  for  u 

and  X  .  The  non-analytic  techniques  include  (l)  difference  approximations 

(2)  Lagrangian  interpolation. 

3.1  Introduction 

Having  established  that  the  moving  boundary  problem  can  be 
expressed  in  the  form  of  non-linear  integral  equations  and  having  realized 
that,  in  general,  these  integral  equations  are  not  expressible  in  terms  of 
known  functions,  we  seek  numerical  solutions.  We  consider  Neumann's  solution 
which  assumes  that  X(t)  a  t^2,  X(t)  being  the  boundary.  Evans  II,  Isaacson, 
and  MacDonald  (1950)  give  three  solutions,  the  first  of  which  employs  Laplace 
transforms,  the  second  assuming  a  power  series  for  X(t)  and  the  third 
assuming  a  power  series  for  u(x,t)  .  Boley  (I96I)  also  uses  analytic 
techniques  to  obtain  solutions  for  a  few  problems. 

Included  in  the  non-analytic  techniques  are  difference  equations. 
Douglas  and  Gallie  (1955)  give  a  system  of  A.  E.  for  which  they  prove 
convergence  and  stability.  Trench  (1959)  obtains  similar  results.  Crank 
(1950)  considers  a  A.  E.  scheme  as  well  as  a  method  using  Lagrangian 
interpolation. 

We  will  examine  the  techniques  and  indicate  some  instances  of 


their  applicability 


.V- 3.f(  ii-.: 
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(3.2)  Analytic  techniques 


(3.2.1) 

Neumann's  solution 

The  assumption  which  Neumann  makes  is  that  the  equation  for 

the  moving  boundary  is  X(t)  a  ,  We  illustrate  this  technique 

for  several  problems.  Each  problem  satisfies 


(3. 2. 1.1) 

=  Ug(X(t),t)  »  Tx 

(3» 2. 1.2) 

Kl“lx  ‘  K2U2k  ■  WXC 

(3. 2. 1.3) 

kiuixx  -  un  x  s  x(0 

(3. 2. 1.4) 

k2U2xx  =  V  *  ^  X(t) 

where  K^,  k^,  and  p  are  thermal  constants.  In  addition  to 
the  above  conditions,  there  are  conditions  which  must  be  satisfied  at 
fixed  boundaries. 


Problem  1 

Melting  for  x  <  0  . 

(3.2. 1.5) 

Ug  -»  V  as  x  -»  oo  f 

(3.2. 1.6) 

u^  =  0,  x  =  0  . 

3£ 

Letting  y  =  -pr  ,  a  particular  solution  which  satisfies 
s/t 


(3.2. 1-3) 

and  (3. 2. 1.6)  is 

(3.2. 1.7) 

u.  =  A  erf  f  — ^ —  J 

1 

•■'O.i  :‘;i> ~a'  ;  ■■ 
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A,  a  constant.  If  B  is  a  constant, 


(3. 2. 1.8)  u0  =  V  -  B  erfc  ( -JL- ) 

2  Wk2  J 

satisfies  (3. 2. 1.4)  and  (3. 2.1. 5)  .  Then  (3. 2.1.1)  requires  that 


(3.2. 1.9)  A  erf 


X(t) 


V  -  B  erfc 


m  1  .T 

-  2(k2t)l/2J  1 


Since  (3. 2.1. 9)  must  be  satisfied  for  all  values  of  t,  X  must  be 
proportional  to  t  '  2  ,  that  is 


(3.2.1.10)  X(t)  =  2*  (kxt)V2 


where  A  is  a  constant  determined  by  (3. 2. 1.2)  .  Using  (3*2. 1.7), 
(3. 2. 1.8)  and  (3. 2.1.10)  in  (3. 2. 1.2)  we  have 

(3.2.1.11)  A  exp  (“X2)  -  KgB^/kg)*/2  exp  (-kjA2/kg)  =  7<Lk^pJii 

We  use  (3. 2. 1.9)  and  (3*2.1.10)  to  eliminate  A  and  B  from  (3*2.1, 
and  solve  for  A  .  Then 

erf  (  -XJ) 

Wk  / 


ll  1  erf  (X) 


(V-T. )  erfc  ( 

1  Wk, 


u2  =  V- 


erfc 


k  J/2 

Hi 


The  solution  of  (3*2.1.11)  for  A  requires  the  numerical  evaluation 


11) 


of  erf  . 


-  8:  -  - 
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Problem  2  Supercooled  liquid 


Suppose  that  the  melting  point  of  a  solid  is  Tp  that  the 


region  x  >  0  initially  contains  liquid  at  temperature  V  <  Tp  and 
that  solidification  starts  at  the  plane  x  =  0  and  moves  to  the  right, 
no  heat  being  removed  from  the  solidified  material  whose  temperature 
will  thus  have  the  constant  value  throughout.  This  problem  takes 

the  following  form: 


(3.2.1.12)  u^(0,t)  =  Tx  >  u2(0,t), 


(3.2.1.13)  Ul(x,t)  =  T1#  x<X(t) 


(3.2.1.14)  u  (x,0)  =  V,  x  >  X(t)  . 


A  particular  solution  satisfying  the  above  conditions  is 


(3.2.1.15)  ut  =  T  ,  x  <  X(t) 


x  >  X(t)  , 


A  being  a  constant.  If 
x  =  x(t )  give 


then  the  conditions  at 


V  +  A  erfc  (A)  = 

A  exp  (~A2)  =  Al  rt1/2  . 


Then  A  is  the  only  zero  of 


S  O  if?  .  '■  .  Hr:  r  >-  ;  j 

,  •  ■  ■  ,  1  ii  '  tr. 

s**Jq  ,  n  ,  a  bite  ill  si  ufc 


'  '  -■■■■'  ■  '  '  :■  so  .  v.  ;  Uv 

.  .  ••:  :■  .■  :  :  V' 


t  '  :;v 


’  (s)jI  < x  .  '  ••  u  \  .  : 
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2  (T1  ‘  V) 

A  exp  (A  )  erfc  (A)  =  - - z-t?  • 

L  * f 


Problem  3  Melting  for  x  >  0 


(3.2.1.17)  u2(x,0)  =  0  , 


x  >  X(t) 


(3.2.1.18)  Ul(0,t)  =  u2(0,t)  =  V 

If  X(t)  =  2A(k2t)1/2  ,  we  have 


A  being  determined  by  the  initial  and  boundary  conditions.  This  sol- 
ution  is  the  same  as  that  for  problem  1  with  the  role  of  the  thermal 
constants  of  the  two  phases  and  those  of  and  V  -  T^  interchanged. 

The  above  problems  should  illustrate  the  method.  These  are 
the  most  important  exact  solutions  where  the  region  x  >  0  is  initi¬ 
ally  at  a  constant  temperature  V,  greater  than  the  melting  point. 
However,  there  are  no  closed  solutions  for  other  important  boundary 
conditions  at  x  =  0  such  as  constant  flux  or  radiation.  The  exact 
solutions  are  frequently  useful  as  starting  solutions  for  other  problems. 
In  such  cases  we  resort  to  other  techniques  to  obtain  solutions. 


r:.,..- ' :  :■  •-."ors  uix..-  ,  ~;>Jt  7  /  ;  ;  IS  ;s  »--u\e.  x.xi  -■  - 

'■  '-b  -■  • .. , : a  .  .9.  ■  v.  >•;  . 
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3.2.2  Evans  II,  Isaacson  and  MacDonald  on  Stefan  problems 

Evans  et  al.  (1950)  give  three  methods  of  solution  (l)  Laplace 
transform  (2)  power  series  for  X(t)  and  (3)  power  series  for 
X(t)  assuming  u(x,t)  has  a  power  series  expansion. 

We  consider  the  problem 

(3.2.2. 1)  kuxx(x,t)  =  ufc(x,t)  0  <  x  <  X(t) 

(3. 2. 2. 2)  u  =  0  x  >  X(t) 

with  boundary  conditions 

(3. 2. 2. 3)  AX.(t)  -  u  (x(t),t)  ,  A,  a  constant 

(3. 2. 2. 4)  X(0)  =  0 

(3.2.2. 5)  ux(0,t)  =  g(t)  . 

We  relate  the  solution  obtained  by  the  use  of  F  and  G  functions  of 
the  formal  solution  with  a  Laplace  transform  solution. 


The  condition  at  x  =  0  is  satisfied  by 


p*  h(t1)  /  _x2  \ 

Jo  1^7^  exp  fa*?) 

/  h(t^)  E(x,tJ  Q,t^)  dt^ 

J  r\ 


dt. 


(3.2.2 .6) 


■i.., 


•  :  :  S’  ’  .  :i3i  ■,  :  ■■  v  •.  ; 


- 
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k\V2 


where  h(t^)  =  ^  ^  J  g(t^)  and  this  satisfies  (3.2.2. l)  .  To  this 

we  add  functions  F  which  do  not  disturb  the  boundary  conditions 


(3. 2. 2. 7)  u(x,t)  =  f  f  (tx)  (E(x,kt;  X(t1),kt1)  +  E(x,ktj  -Xfcj),^))  dt][ 

J  n 


/  1,  \^/2  nt 

-  ^  -  J  J s(tx)  E(x,t;  0,tx)  dtx  . 


Differentiating  with  respect  to  x,  we  have 


ux(X(t)  -  0,t)  =  AXt(t) 


„  1/2  ,  rt,  £(t  ) 

'  -  S  f 

J  o  1 


E(x(t),  kt;  x(t1),kt1)  +  (x(t)  +  x(tl)) 
E(x(t),kt;  -  x(t1),kt1)3  dtx 


2(x/k)1/2  fo 


t  g(tx)  x(t)  E(x(t),t;  o,ti; 


dt- 


We  can  extend  this  solution  to  the  right  of  X(t)  .  Since  u  =  0 
for  x  >  X(t)  , 


(3. 2. 2. 9) 


ux(x(t)  +  0,t)  =  0 


and 


(3.2.2.10) 


0  =  -(*/k)1/2  £(t)  -  i 


o  2(<t/k)1/2  Jo 


. 
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Subtracting  (5. 2. 2. 8)  and  (3.2.2.10)  we  have 

(k/*)1/2  AX  (t) 

f(t)  .  - 

and  u(x(t),t)  =  0  give  the  following  integral  equation  for  X(t); 

(3.2.2.11)  2  j  g(t  )  E(x(t),t;  0,^)  dt1 

J  o 

=  A  f  Xt(t)(B(X(t),kt;  x(t1Jjtt1)  +  E(x(t),kt;  -X(t  ),kt  ))  dt 
J  o 

This  solution  satisfies  (3. 2. 2.1)  for  x  <  X(t)  .  For  x  >  X(t)  ,  u  =  0 
along  a  curve  X^(t)  arbitrarily  close  to  X(t)  since  u  is  continuous 
on  crossing  X(t)  and  we  have  made  u  =0  along  X.(t)  .  Thus  u  =  0 
for  x  >  X(t)  . 

Evans  et  al.  obtained  this  solution  in  another  way;  they  used 
Laplace  transforms.  Inverting  x  =  X(t),  we  define  the  boundary  by 

t  =  f(x)  . 

The  region  x  >  X(t)  is  defined  by  t  <  f(x)  since  the  curve  passes 
through  the  origin  and  in  this  region  u  =  0  .  The  Laplace  transform 
u  is  defined  by 


(3,2.2.12) 


e"pt 


udt  . 


and  u 


xx 


do  not  exist  on  the  boundary  but  we  can  define  the  limits 


i(  0  *:)s  -  (  ..}>■'  :-y.i  '  .  • .  •  )./ 

. 

■-'■  ■■  '  :  '  ,  -i.'.  , 

E 


■  -  'V-  oififs  r:  .  h;r-  tij  ^rjj  j 

beoilab  aj 


v  ■•••  *•' ■  !■  '  ai;<-'  vo. .  • .  i 
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e  ^>t  u  dt  and  /  e  u  dt  e  >  0 

f(x)+e  1  Jf(x)+<=  xx 

In  this  way  the  D.  E.  for  u  is  transformed  to 
Uxx*pu/k  =  -Ae 

and  the  boundary  conditions  become 

uxx(°>p)  =  S(p)  ,  u(oo,p)  =  0  . 

The  solution  of  this  equation  followed  by  an  inversion  from  the  parameter 
p  to  t  gives  us  the  integral  equation  (3.2.2.11)  . 

The  above  is  an  ingenious  use  of  the  Laplace  transform.  However, 
the  use  of  Laplace  transforms  assumes  the  definition  of  u  and  its  partial 
derivatives  for  all  t,  0  <  t  <  °o  and  the  convergence  of  the  integrals 
(3,2.2.12)  and  (3. 2. 2.I3). 

A  power  series  solution  for  x  =  X(t)  of  (3.2.2.11)  may  be 
found.  Let 

CO 

(3.2.2.14)  x(t)  =  c.t1 

i=l 

and  assume  g(t)  =  g  ,  a  constant.  c^  and  c^  are  determined  using 

(3.2.2.14)  to  expand  integrals  of  (3.2,2.11)  in  powers  of  t  .  In  the 
limit  as  t  tends  to  zero  (3.2.2.11)  becomes 

4> 

cL  =  g/A  . 

2 

If  we  replace  X(t)  by  gt/A  +  c^t  and  X(t)  by 


• 

■  :  ■  '  -.S;  ;;  i. :  .  •.■■■*  .  •  •  .  w.,  ■ 

'  ■-!•£  ■  iT.>  |  .  ,  •  v  ■  • 

■  ''  &  ■■ 

3-  '  N  iJ'  S.i, c  .:;.br :  .  •;  •  z  ■  \.i 

"i ;  ■ n  •  *  y  J  '  <  .  -i' 
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t  X  (0) 

x(t)  =  x(o)  -  txt(o)  +  — y\ — 


then  (3.2.2.11)  becomes 


4gt 


1/2  2g2  n1/2  t  g3 

aA  2  2 

a  A 


1/2  16Ac  t3//2 
=  bgtL^  +  - 1 - 


2g2  It1/2  t  ng5  t5/2 

aA  +  2  2 

3a  A 


giving 


2a2A3 


The  determination  of  further  coefficients  by  this  method  becomes  extremely 
tedious . 

Evans  and  his  co-authors  obtained  another  power  series  repre¬ 
sentation  of  X(t)  by  assuming  u(x,t)  has  power  series  about  (0,0) 
for  all  u  in  0  <  x  <  X(t)  .  They  assume  that  any  discontinuities  in 
the  derivatives  of  u  occur  only  on  crossing  X(t)  . 

By  differentiating  (3. 2. 2. 2)  and  (3.2.2. 3)  along  x  =  X(t) 
and  of  (3.2.2. l)  in  0  <  x  <  X(t)  and  using  (3.2.2. 5),  all  derivatives 
of  X(t)  may  be  evaluated  at  the  origin.  If  in  (3. 2. 2. 5),  g(t)  =  g  , 
constant,  then 


(3.2.2.15) 


x(0«S£ 


g' 


2!kA' 


5r 


3!kA; 


,3 


51g7  >  x  827g9  .5 

v  3  7  C  +  4  9  C 

4! k  A  5 ! k  Ay 


Further  coefficients  of  powers  of  t  may  be  obtained  using  the  following 


bo  if  7-305'  elite  •/.-?  ;  .  na.  rlii  -  ro  ^  .  >rl' 


a  si  r a  tewoq  a  7  "I  j,  .  .o  g:olr’.:  vJ  )•' 

<  ,*X  sol, a  07  0  v.  •  J"3r>o 


- svi 3 K  i/7  '  3 b  I  r B  <  ( ,S,S  '  g rr .  b •  B 
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relations 


we  have 


and 


where  a^ 


g  =  g(0, 


where 


Letting 


V  i  j 

u  =  )  a. .x  t 

L 


00 


i,j=0 


and 


X(t)  =  y 

/  i 

i=0 


V  ' 


00 


X  <1+1>  ai+llJ  (  I  v" 

i,j=0  n=0 


00 


=  A 


y  (n+l) 
n=0 


c  ,  t 
n+l 


00  00 

I  ■»(  I  V ) 1,1 

i,j=0  n=0 


=  0  , 


a . 


i+2,  j 


(i+l) (i+2)a 


S  ai >  j+1 


.0  =  8  and  c0  =  0  • 


If  g  is  no  longer  a  constant  but  an  analytic  function  of  time, 
X(t)  is  found  to  be 


°o  i 

c .  t 


*•>  ■  i  s 

iTl 


C1  = 


g(0) 

A 

 g3(o) 


2.3' 

a  A 


gt(o) 


A 


...  '  •  ■  ■  "  ■  • 


:  .  . 


■ 


..  '  '  i  '  il" 
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C  -  ?g?(0) 

5  ”  2.5 


6g  (o)g2(o) 


.V 


*t<°) 

A 


and  so  on. 


Evans  et  alii  consider  the  case  where  the  material  of  uniform 
thickness  has  infinite  length.  The  problem  is  given  by 


U2t  a2U2xx 


0  <  x  <  X(t) 


u2(X(t),t)  -  0  .  ux(X(t),t) 


u..  =  a.u, 

It  1  Ixx 


X(t )  <  x  <  oo 


with  boundary  conditions 


Kgu2x(X(t),t)  -  KlUlx(X(t),t)  =  pHXjt)  =  BXt(t)  , 
u^(x,0)  =  f(x)  where  f(x)  <  0  and  f(0)  =  0 
*2x(0,t)  =  g(t)  , 

where  p  is  the  density,  H  the  latent  heat,  and  the  coefficients 

of  thermal  conductivity. 

Tn  this  problem  we  assume  g  is  analytic  in  t  and  f  has 
a  Taylor  expansion  with  an  infinite  radius  of  convergence  about  x  =  0  . 
Using  the  method  of  solution  of  the  previous  problem,  we  have 


x(t)  = 


K2g(0)  -  Klfx(0)  fc 


B_ 

2! 


-K  m  (Kas(°)  -  Kifk<° 
1  l*2  B 


■ 


:  '.r 


-  '•  1  ‘V:i'  '^:;r 


■■  if. 


:  >i: ,  "v;  m;:  ;  f  \>  a  •  Ml  cl 


-  38  - 


Upon  examination  of  the  series  for  X(t),  we  can  see  the  degree 


of  agreement  of  the  first  five  terms  of  X(t)  in  (3.2.2.15)  .  However, 
we  note  that  the  method  of  determining  coefficients  of  powers  of  t 
in  X(t)  becomes  extremely  tedious. 

3.2.3  Boley's  method 

Boley  (I96I)  transforms  the  partial  differential  equation  to 
integro-dif ferential  equations  as  did  Evans  and  his  co-authors.  To  solve 
these  equations  he  finds  an  analytical  solution  valid  for  short  times 
and  then  proceeds  numerically. 


We  consider  one  of  the  problems  to  illustrate  the  method. 


Kuix(°,t)  =  -qQ 


0  <  t  <  t 


m 


This  problem  yields  the  integro-dif ferential  equations 


m 


■  J  ... 


I  ■  -  ■  '  ■  '  •  • 
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(3.2. 3- 2) 


t-t 

■X(t)  f 
J  o 


m  f(t-t..) 

— -  E(x(t),kt  ;  0,0)  dt 

1 


2(k7t)1/2  -{Qq  erf 


{Q° erf  (w^)  '  6  LX‘(t) } 


under  the  condition  X(tm)  =  0,  where  f  is  unknown,  is  the  melting 
temperature  and  t^  is  the  time  the  surface  reaches  melting  temperature. 
Changing  to  the  following  dimensionless  variables 


t  Q0x(t) 

=  —  -  1  ,  XL(t)  =  ~l/2 

m  it  <  l 


«V£cT 


it' '  kT. 


’  m  “  2L 


m 


m 


F(y)  =  ~fQ^'  > 

w0 


(3. 2. 3,1)  and  (3. 2. 3. 2)  become 


ry 

(3.2.30)  /  F(y-yx)  E(X1(t),y1;  0,0)  dyx  = 

j 


2  ((1-(lt(Uy))V2)ier£c  )) 


(3.2. 3.U)  X 


.«  L 


y  F(y-y^)  E(xi(y),y1;  0,0) 


dy. 


=  *V2 


erf 


(Xl(y))  2Xly 


(i+y) 


1/2 


m 


The  condition  X(tm)  =  0  becomes  X^(0)  =  0  .  It  follows  that 


'  •  .v  ..  n  ..  .  1  ■  ; 


stt 


■  :  ,i. 


(  ( 


. 


■ 
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(3.20.5)  lim  (  ~7i')  =  lim  r-Hfcl)  =  0  for  n  <  1  , 

and  that  E(X(y),y^;  0,0)  for  0  <  <  y  is  very  close  to  1  at  small 

values  of  y  for  nearly  the  entire  range  of  integration,  and  drops  off  to 
zero  near  the  lower  limit.  The  right  hand  side  of  (3-2. 3.3)  expanded  for 
small  y  is 

-y  +  +  2*1/2  Xx(y)  . 

If  y  is  sufficiently  small,  only  the  first  term  need  be  kept  in  view 
of  (3, 2. 3.5).  Thus  (3* 2. 3, 3)  becomes 

py  FCy-yJ 

(3.20.6)  /  — jj2~  dyi  =  -  y  y  <  <  1  . 

Jo  yx  ' 

This  is  a  special  case  of  Abel's  equation  which  has  the  solution 

.  1/2 

(J.2.3.7)  F(y)  =  . 

Substituting  into  the  left  hand  side  of  (3.20.4)  we  get 

2Xi  py  (y-yJ1/2 

J  E(X1(y),yi;  0,0)  - dy± 

=  -2y1/2  ierfc  (X^y1/2) 


-2(y/x)^2  +  2Xj  + 


41 


The  right  hand  side  of  (5,2, 3,4)  with  (3,2, 3, 5)  for  small  y  is 


2k 


1/2 


X,  + 


m  ly 


Thus 


(3. 2. 3.9) 


1/2 

X,  =2Z_i. 
ly  k 


or 


V*>  -  H 


3/2 


3c 


Equations  (3.2. 3*7)  and  (3*2. 3. 9)  represent  a  short-time  solution^  to 
extend  its  range  it  is  now  assumed 


F(y) 


-gy1/2 

m 


Z  V 

i=l 


r(i+l)/2 


and 


xl(y)  -  SH 


3/2 


OO 


3c 


z  V 

i»2 


(i+2)/2 


The  coefficients  are  obtained  upon  substitution  into  (3*2. 3*3)  and  (3.2. 3,4) 
and  equating  coefficients  of  like  terms. 

An  estimate  of  accuracy  may  be  obtained  by  comparing  the 
contribution  of  the  last  term  considered  with  the  sum  of  all  the  terms. 

Boley  then  obtains  a  numerical  solution  for  y  >  y*  with 
sufficient  accuracy.  For  y  <  y* 


fy 

/  F(y-yj )  E(X1(y),y];  0,0)  dy^ 
Jr* 
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y-y* 

/  Ffr-y^  E(X1(y);y1i  0,0)  dyx 

J  o 

-  |  J  (y/y^i)1/2  exp  {-\/yx)  dvx  «  +  i2  . 

y-y* 

I2  can  be  evaluated  readily;  1^  is  approximated  using  numerical 
integration. 

As  in  Evans*  method,  the  evaluation  of  the  coefficients  in  the 
two  series  is  tedious. 

3.2.4  Summary 

It  is  apparent  that  analytic  techniques  are  unsatisfactory 
because  they  are  tedious  ~~ repetitively  evaluating  coefficients  of  a  power 
series  by  equating  like  terms  in  an  equation  --  and  no  radious  of  con¬ 
vergence  has  been  given.  We  would  like  solution  techniques  for  machine 
computation  rather  than  semi-analytic  techniques  whose  precision  is 
undefined . 

3.3  Non-analytic  techniques 

Having  noted  that  the  analytic  procedures  in  the  previous 
section  are  not  machine  type  procedures,  we  seek  non-analytic  procedures 
which  are  machine  type,  Thus  far  the  methods  are  A.  E,  approximations 


to  the  D.  E.  . 


8f,r-T  ol:i  -;x 

. 
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3.3*1  Implicit  difference  system  of  Douglas  and  Gallie 

Douglas  and  Gallie  (1955)  treat  the  problem  defined  by  (2.1.1.1) 

-  (2.1. 1.5)  .  One  can  see  that  (2. 1.1.4)  is  equivalent  to 

rX(t) 

X(t)  =  t  /  u(x,t)  dx 

J  o 

In  their  paper,  they  introduce  an  implicit  A.  E.  system  equivalent  to 
(2. 1.1.1)  -  (2. 1.1.5)  and  demonstrate  that  its  solution  converges 
uniformly  to  that  of  (2. 1.1.1)  -  (2. 1.1. 5)  in  a  suitably  restricted 
domain.  They  show  their  method  to  be  stable  with  respect  to  small  errors. 
By  stable  is  meant  that  property  P  of  a  solution  which  is  preserved 
when  the  initial  specification  of  the  problem  is  subjected  to  variations 
v  .  In  our  case  the  property  P  is  boundedness  and  the  variations  v 
are  round-off  errors.  We  note  in  passing  that  an  iterative  process 
carried  out  exactly  may  yield  a  sequence  of  functions  which  converge 
uniformly  to  a  desired  function  yet  boundedness  may  not  be  preserved 
when  the  iteration  embodies  round-off  errors.  Their  proof  is  based  on 
a  maximum  principle  and  properties  of  equicontinuous  families  of  functions. 
We  outline  the  procedure  below. 

In  this  procedure  we  fix  Ax,  the  space  step  and  vary  At, 
the  time  step  so  that  the  moving  boundary  is  always  kept  at  the  right¬ 
most  mesh  point.  A  sequence  (At(k)),  k  =  0,  1,  ...  is  obtained 
inductively  and  iteratively.  The  following  notation  is  useds 
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x(i)  =  lAx  , 


n~l 


u(i,j)  =  u(x(i),t(j))  . 


We  let 


u(0,0)  =  0  ,  u(Q,l)  =  A  x  , 

u(l,l)  =  0  ,  At(0)  =  Ax 


and  take 


At(n,0)  >  0  , 

arbitrarily.  We  use  the  following  algorithm  which  is  the  difference 
analogue  of  (2. 1.1,1)  -  (2. 1,1. 3)  where  the  boundary  is  assumed  at 

((n+l)  t 


Au(i,n+l,r)  = 


u(Q,n+l,r)  -  u(l,n+l,r)  =  Ax 
u (n+l, n+l, r)  =  0 


x  «  ®,n 

r  >  0 


t 


At(n,r+l)  = 


u(i,n+l,r))Ax-t (n) 


Thus  (u(i,n+l,r))  and  At(n,r+l)  are  computed  and  the  iteration  is 
continued  to  an  arbitrary  point,  say  q  times,  to  obtain  {u(i,n+l,q) } 


» 
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and  At(n,q+l)  .  Then  we  set 

u(i,n+l)  *  u(i,n+l,q)  , 
and  At(n)  =  At(n,q+l) 

and  proceed  to  the  next  time  level. 

It  is  known  that  the  solution  of  a  parabolic  D.  E.  satisfies 

a  maximum  principle.  However,  this  is  not  true  for  all  A.  E, 

2 

(e.g.  the  explicit  A.  E.  ut  n ^  =  ut  n  +  AtAxUi  n  with  an  un®table 

2 

ratio  At  to  (Ax)  )  .  This  principle  does  carry  over  in  the  above 
algorithm  as  is  shown  in  Douglas  and  Gallie  (1955) •  They  also  show  that 
the  iteration  gives  a  sequence  of  functions  tR(x)  which  converges  to 
t(x)  .  Corresponding  to  each  t^(x)  *-s  X^(t),  the  inverse  of  t^(x)  . 

The  authors  show  that  X^(t)  converges  to  X(t)  ,  the  inverse  of  t(x)  . 
Having  established  the  convergence  of  t^(x)  and  X^(t),  they  show 
that  as  Ax  -» 0  u(x,t,r)  converges  uniformly  to  u(x,t)  . 


Besides  showing  that  their  method  converges,  they  show  it  is 
stable.  For  a  method  to  be  useful  for  computation  it  must  be  stable 
against  round-off  or  other  small  errors.  Suppose  we  perturb  t(n),  At(m-l) 
and  u(i,n)  along  some  line  t  =  t^,  n  >  p  .  Demoting  the  perturbed 
values  by  t(n,*)  ,  At(n-1,*)  and  u(i,n,*)  we  have 


d(i,n)  =  u(i,n,*)  -  u(i,n) 


where  max  d(i,n)  <  d  .  Then 
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A^d(i,n+l)  = 

X 


u(ijrn+lj*)  -  u(i,n,*) 

At(n,*) 


At(n) 


dl1i,n+1),r-.  d(iini 

At(n>) 


.  «{i*s±LL . -  “J&iil  /MslL, 

At  (n)  \At  (n,*) 


'  At(n) 

Atlt;*) 


or  d(i,n+l)  =  d(i,n)  -  A  u ( i, n+1 ) (At (n)  -  At(n,&))  + 

<o» 


At(n,*)  Axd(i,n+!) 


n-1 


Now  by  t(n)  -  t(n,*)  =  - d(i,n)Ax 

i=l 


we 


obtain  At (n)  -  At (n?*)  =  d(i,m+l)  Ax  ~  (t-t(n;*))  . 

i=l 


By  extending  their  proof  of  their  maximum  principle  they  conclude 

max  jd(i,n-s-l)  j  <  max  jd(i,n)  j  +  jAt(n)  -  At(n#*)  | 

i 

Now  j t: ( n)  -  t(n,*)  j  <  (n~l)Ax  max  |d(i,n)  j 


and 


|At(n)  -  At(n,*)  j  <  nAx  max  jd(i,n)  j  +  |t(n)  -  t(nJ1*)  j 


max  jd(i,n)  | 
i 


<  lHr(n-"l)Z^X 
—  1-nAx 


max  |d (i, m) j 

i 


Thus 
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Hence, 


H“1 

max  |d(i,n)  |  <  d  n 

k=p 


l-i-(k-l)Ax 

1-kAx 


Now  for 


and 


x  <  1  -  6  ,  n  <  (l-e)/Ax, 

2  (l~e)/Ax 
max  |d(i,n)|  <  (-) 

| t (n)  -  t(n,*) |  <  D  . 


d 


ss  0 


Since  Ax  is  considered  fixed  in  this  argument,  no  error  introduced  into 
the  calculation  can  grow  unboundedly.  Thus,  the  method  is  stable. 


The  authors  made  calculations  for  Ax  *  0,1,  0,05,  0,025 
using  16  iterations.  They  noted  that  t  (x)  increased  as  Ax 
decreased  and  the  convergence  was  apparently  linear.  Upon  using 

At (n,0)  =  2At(n»l,0)  -  At (n-2,0)  , 

an  additional  iteration  yielded  4  significant  figures  in  At(n)  and 
u(i,n+l)  and  5  significant  figures  in  t(n+l)  , 

It  would  seem  that  fixing  q  ,  the  number  of  iterations  is 
an  unnecessary  restriction  on  the  problem.  It  appears  that  one  could 
specify  a  tolerance  which  would  be  acceptable  and  proceed  iteratively 
until  the  values  of  u  are  obtained  to  the  desired  degree  of  tolerance. 


5.3*2  Results  of  Trench 

Trench  (1959)  formulates  an  explicit  numerical  analogue  which 
is  a  finite  difference  technique.  He  shows  that  the  mesh  functions  so 
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defined  converge  to  X(t)  and  u(x, t)  and  that  convergence  is  uniform 
in  every  finite  interval.  He  shows  that  existence  is  a  consequence  of  the 
mesh  functions  and  an  existence  theorem  for  a  classical  boundary  value. 
Trench's  convergence  and  stability  results  are  essentially  those  of 
Douglas  and  Gallic  (1955)*  However,  Trench  gives  an  explicit  method 
whereas  Douglas  and  Gal lie  give  an  implicit  method. 


The  problem  considered  is 


=  ut(x,t) 


ux(0,t)  =  f(t) 


u(x(t),t)  =  g(t) 


where  f  and  g  are  such  that 


dx 


for  which  he  gives  the  algorithm 


m  =  1,  2,  .  . , ,  M  “1 


(3. 3. 2.1)  ^ 


2Xu(Mn-l,n) 


u  (M  ,  n-j-1 )  = 
n 


V 


+  KO(h3) 
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Also  if  Mn+i  33  *-+Mn  which  occurs  if 


X(n)  <  (Mn-}5)h  <  X(n+1) 


then 


/„  ,  x  u(Mni’Tl+l)  g(n+l)  ,rrs/t  2S 

u(m  n+i)  =  p  T7p~T~  +  up“?  +  K°(h  > 

n+1  n+1 


Now, 


M  = 
n 


W  -  *  ■ 


X(n)  -  M  h 
'  n 


n 


where  h  and  k  are  the  increments  in  x  and  t  respectively  and 

p 

1  >  h  >  0  ,  K  =  Ah  ,  A  a  constant  such  that  0  <  A  <  1/2  „ 

That  portion  of  (3.3»2.l)  which  is  homogeneous  in  u(m,n)  can 
be  formulated  more  conveniently  as  follows? 


S  is  a  (l+M  ) “dimensional  vector 


S  = 


s  (0) 

s(l) 

s(2) 


s(Mn) 


and  H(n)  is  the  operator  on  S  which  produces  the  (l-iM  ) -dimensional 


1 


[x]  means  the  greatest  integer  less  than  or  equal  to  x  . 
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vector  R 


H(n)S  with 
r (O)  -  r (l) 

r (m)  =  As  (m-1 )  +  (l-2A)s(m)  As(m+l), 


2As (Mn-1) 


m  ^  la2a  teej  “  1 

n 


and  if 


M  .  =  1:-M 

n-fl  n 


r -Mn+1'' 


P  ,r(M  ) 
n+1  '  n' 


1+P 


n+1 


H(n)  is  a  linear  operator  dependent  on  the  sequence  (X(n) ]  and  the  mesh 
size  h  . 

Trench  shows  H  to  be  a  bounded  operator,  its  norm  bounded  by  that 
of  S  .  In  his  paper  Trench  also  proves  unicity  of  the  solution. 

3. 3*3  Two  methods  of  Crank 

Crank  (195?)  gives  two  methods  for  numerical  treatment  of  heat 
flow  problems  in  which  a  transformation  boundary  moves  through  the  medium. 

The  problem  is  transformed  by  a  change  of  variable  from  one  involving  the 
simple  heat  equation  with  an  awkward  moving  boundary  to  an  eigenvalue 
problem  in  which  the  equation  is  slightly  more,  complicated  but  the  bound¬ 
ary  is  fixed.  In  the  second  method  Lagranglan  interpolation  formulae  are 
used  to  develop  finite  difference  approximations  to  space  derivatives  of 
temperature  based  on  values  at  points  unequally  spaced  in  the  direction 
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of  heat  flow.  This  enables  one  to  follow  the  course  of  the  transformation 
boundary.  We  proceed  to  examine  the  two  methods. 

To  illustrate  the  method  of  transformation  of  variables  Crank 
considers  the  diffusion  of  solute  in  solution.  The  corresponding  problem 
of  heat  flow  in  a  medium  is  that  in  which  latent  heat  of  transformation 
is  absorbed  or  liberated  at  a  definite  temperature.  In  the  diffusion  we 
consider  an  infinite  slab  of  uniform  material  of  thickness  2a  placed 
symmetrically  in  a  solution  of  extent  2L  .  Solute  is  allowed  to  diffuse 
into  the  slab  which  is  initially  free  of  solute.  The  slab  has  a  number 
of  fixed  sites  on  each  of  which  one  and  only  one  diffusing  molecule  can 
be  instantaneously  and  irreversibly  immobilized  and  so  withdrawn  from  the 
diffusion  process.  There  are  S  sites  per  unit  volume  of  the  slab;  the 
concentration  of  solute  in  solution  is  always  uniform  and  is  u^  initially, 
expressed  as  the  number  of  molecules  of  solute  per  unit  volume.  Denoting 
by  u  the  concentration  of  freely  diffusing  molecules  at  any  point  in  the 
slab  in  number  of  molecules  per  unit  volume,  the  concentration  of  molecules 
immobilized  on  the  sites  is  zero  if  U  =  0  and  is  S  when  U  >  0  , 


Taking  x  to  be  the  distance  coordinate  measured  in  the  direction 
of  diffusion  from  X  =  0  at  one  surface  of  the  slab  and  using  the  dimen¬ 
sionless  variables 


x  =  X/a,  t  =  DT/a  ,  u  =  U/UL  ,  s  =  S/U 


where  D  is  the  diffusion  coefficient,  we  have  to  solve 


ufc(x,t)  *  uxx(x,t) 


(50.3.1) 


XT  U 


0  <  x  <  x(t) 


1 

,  , 


- 


i 
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(3. 3. 3.2)  u(x,0)  =  0 


0  <  x  <  1 


(3. 3. 3. 3)  u(x,t)  =  0 


x  >  X(t ) 


*'  o 


t  >  0 


(3. 3. 3.6)  X(0)  =  0  . 

X(t)  is  the  value  of  x  at  which  u  =  0,  and  for  all  x  <  X  all  sites 
are  occupied  and  for  x  >  X  they  are  all  unoccupied.  X(t)  increases 
from  0  as  the  boundary  between  occupied  and  free  sites  moves  into  the 
slab . 


The  singularity  on  the  surface  of  the  slab  x  =  0  is  handled 


by  introducing  new  variables  p  =  x/t ^/2  ,  q  =  t1/2,  and  (3*3»3«l) 
-  (3.3o3.^)  become 


(3. 3. 3.7) 


(3. 3. 3. 8)  u(»,0)  =  0 


(3.3. 3.9)  u(p,q)  =  0 


P  >  Px  >  q  > 


where  P  =  X/t1/2  =  X/q  and 


(5.3.3.10)  -P  /2  =  up( P,  t )/ sq-i  P/2q 


while  (3. 3* 3«5)  becomes 
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(3.3.3.11) 


qdp 


(3.3.3.12)  2u  (p,0)  =  PU  (p,0) 

r  r  tr 

(3.3.3.13)  u(0,0)  =  1 

(3*3.3.14)  u(p,0)  =  0 

(3.3.3.15)  2up(P#q)  =  -sP  . 


In  deriving  (3.3*3.15)  from  (3«3«3.10),  P  (0)  is  assumed  finite*  A 
solution  of  (3.3.3.12)  subject  to  (3.3*3.13)  is 


(3.3.3.16)  u  =  1  +  B  erf  (p/2) 

where  B  is  a  constant  which  together  with  p  is  determined  by  (3.3*3.14) 
and  (3.3.3.15)  since  we  have 


(3.3.3.17)  1  +  B  erf  (P/2)  =  0 

(3.3.3.18)  Jt1/2  P  +  erf  (P/2)  exp  (P/2)2  =  2/S 


We  fix  the  boundary  by  setting 

(3.3.3.19)  r  =  p/P  ,  q  =  q 
to  pbtain 


(3.3.3.20) 


2u  (p,q) 

i  i 


p2 

p  q« 


+ 


2ur(p,q)  ur« 


>q) 


s 


0  <  r  <  1 


. 

-  i 

■ 

(S\q)  .  tp  S'  4  1  «  u  (fcl ♦£.£«,£) 


,  •  J  ■  , 

,.q 
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The  problem  is  now  that  of  solving  (3.3»5»19)  for  trial  values  of  P  and 
ur (P,q)  which  are  consistent  with  (3»3.3«2Q)  until  (3»3»3«2l)  is 
satisfied . 

Crank  uses  a  finite  difference  technique  to  solve  (3*3. 3*20) 

-  (3»3»3*22)  based  on  the  following  equations  -  - 


=  (p2(m)  +  P2 (m+1 ) )  (u(m,n+l)  -  u(m,n)) 


-  2( (u(m+l,n+l)  +  u(m+l,n))  ~  2( (u(mJ!ln+l)  +  u(m,n)) 


(m-l,n))  , 


A(m,n)  =  ur(m,n) 


-  u 


ur(m+l,q)  =  ^  (nP2(m)  -  P(m)  P(m+l)  -  (ml)  P2(m+l))  -  u  (m,q) 


(3.3.3.24) 


Values  of  u(mJ)0)  are  calculated  using  (3.3. 3. 16)  -  (3.3.3.19)  .  Successive 


steps  in  the  solution  for  u(m,n)  to  u(m,n-;-l)  are  obtained  as  follows: 


u 


■ 

('  -  ),  -  (o  ) 
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(1)  estimate  P(m+l)  at  t  =  (n+l)  At  and  calculate 

ur(m,n+l)  using  . 

(2)  for  m  =  M,  r  =  1,  using 

ur  (m,  n+l )  = 

since  u(M,n)  =  0  with  m  =  M-l  in  (3»3*3»23)  we  have 
two  equations  from  which  to  calculate  u(M-l,m)  and 
u(M-2,n) 

(3)  integrate  backwards  along  t  =  (n+l)At  using  (3-3-3-23) 
to  obtain  u (m-l, n+l)  knowing  u(m,n+l)  and  u(m+l,n+l) 
for  successive  values  of  m 

(k)  adjust  the  value  of  P(m+l)  until  the  value  of  u  at 
r  =  0  agrees  with  that  obtained  using  (3»3» 3*22)  . 

(5)  continue  in  successive  steps  At  until  qP  =  1 

Crank  illustrates  the  value  of  this  method  by  giving  numerical 
results  for  two  steps  in  time  and  four  steps  in  distance. 

The  plan  in  Crank's  second  method  is  to  extend  the  finite 
difference  formulae  used  to  approximate  derivatives  in  such  a  way  that 
they  can  be  used  to  follow  the  motion  of  the  transformation  boundary. 

This  means  developing  difference  approximations  based  on  function  values 
which  are  not  equally  spaced  in  the  argument.  Using  the  Lagrangian 
interpolation  formula 


. 


30  v'., 

.  •  \  .  :  Xy.  ■' 
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u 

f(«)  »  ^  lj(x)  f(Hj) 

j=0 


where 


P„(x) 

lJW 


P^(aj)(x-aj) 


and 


Pn(x)  =  (x-a0)(x-ai)  .....  (x-an_1)(x-an) 


and  restricting  n  to  2  ,  we  find 


(3.30.25) 


f  (x) 


XX 


■  z 


f(ai) 


L  *  (vV 

i=0  j=fi  1  J 


and 


(3.3.3.26)  fx(x)  =  1[  (x)f(a±) 

±A 


where 


400  - 


(x-ai+i)  +  (x-a1+g) 

ai'ai+l^ai'ai+2' 


the  subscripts  being  taken  modulo  3  « 


We  apply  these  formulae  in  the  neighborhood  of  the  moving 
boundary.  We  assume  the  slab  is  composed  of  N  layers  each  of  thickness 
h  and  let  the  boundary  at  T  be  somewhere  in  the  (n+l)st  layer  such 

JL 

that 

X  =  (n-l+p)  h 

where  p  is  fractional  and  1  <  p  <  2  .  If  for  x  <  X  , 


f (a0)  =  u(r-2,t)  ,  f(ax)  =  u(r~l,t)  ,  f(a2)  "  u(x>fc) 


' 
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then  (3»3.3«25)  becomes 


/f(a0) 

u.(r-1't)  =  (— 


f(ax) 


f(a2) 


P  +  P(p+l) 


(Ax)£ 


and  (3*3»3»26)  becomes 


/Pf(a0)  (P+1  f(aJ  (2p+l)  f(a  )  \  x 

u*(x’T)  -  (th— - p -  + 


p(p+i)  "  y  ax 


For  x  >  X  , 


f(a5)  =  u(X,t),  f(a^)  = 


and 


%x(r+1>t) 


f(a3) 

P-2)(p«3) 


u(r+l,t),  f(a5) 


f(\)  f(a5) 

T*  1 


=  u(r+2,t) 


2-p  ‘  3“P 

Combining  these  formulae  and  using  an  approximation  of  the  type 

(3.5.3.27)  HiHUlL- i-uIl^T+AtJ. 

we  get 


A  t 


=  ut(m,T) 


r  ,  .  .  .  2At  /f(an)  ^(ai)  1 (au . 

u(r-l,tHAc)  =  u(r-l,t)  +  - - 5  (-^ - —  +  ^ 


(Ax)£ 


2At  /  T^aj} 


(3.3.3.28) 


u(r+l,t^At)  »  u(r+l, t)  + 


f (a[j_)  f(a5)  ^ 

2-p  +  3-p  ) 


P  -P  = 


2t  /  r(2P-5)  f(a^)  (p-3)  f(a^) 


s(Ax)£ 


-2)  f  (a 


(P“2)(p-3)  P-2 


3“P 


f  pf(ao) 

(p+1)  f(«j.)  (2p+l)  f(a2)  'I  ^ 

t  P+1 

P  ’  p(p+1)  J  ) 

( ~  ■) 
v-  J 


. 

•ffliso’i  tq&  .  3  gj  :  j  bin  . 


+J  *  • 


+  + 


. 
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When  r  increases  by  unity  we  obtain  the  value  of  u(r-l,t)  by  using 


the  Lagrange  interpolation  fomula 


+ 


+  2u (r^jt) 
P  P(P+1) 


and  this  becomes  the  new  u(r-l,t)  used  in  (3*3.3*27)  and  (3*3.3.28)  . 


The  application  of  Lagrange  interpolation  formulae  involves 


many  steps  in  time  but  each  step  is  relatively  simple  to  work.  When  the 
boundaries  are  fixed  by  using  suitable  variables,  far  fewer  steps  in  time 
are  needed  but  each  may  involve  several  trial  solutions. 

3.3»4  Ehrlich's  method 

Ehrlich  (1958)  uses  the  difference  equation  analogue  of  the 
problem  to  obtain  a  solution.  He  uses  a  grid  on  the  space-time  plane  for 
the  following  problem; 


(3.3.4. 1) 


u 


XX 


X(t)  <  x  <  3 


(3. 3. 4. 2) 


0  <  x  <  X(t) 


(3. 3- 4. 3)  u(x(t),t)  =  1 


(3. 3.4.4)  ux(x(t),t)(i~7)  =  xt(t) 


(3. 3. 4. 5) 


u 


(x>0) 


c 


(3.3.4.6a)  ux(0,t)  =  -Pj^f  (t) 


X(t)  =  0 


'.'.J  .  .. 


. 
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(3.3.4.6b)  ux(o,t)  =  -p2f(t)  x(t)  >  o 

(3- 5. 4. 7)  u(0,t)  =  g(t) 

(3. 3. 4. 8)  ux(B,t)  =  0 

where  B  denotes  the  thickness  of  the  material.  Condition  (3. 3*4. 7)  is 
used  only  for  the  non^melting  problem  and  (3.3.4.6a  and  b  )  for  both 
melt  and  non-melt  problems.  With  the  following  configuration 


i-1  i  i+l 

Figure  3.3=4. 1 

and  letting  r  =  ^  ,  the  difference  analogue  of  (3. 3*4.1)  and 

(Ax) 

(3. 3. 4.2)  is 

(3. 3.4.9)  -ru(i-l,j+l)  +  2(l+r)  u(i,j+l)  -  ruCi+l^j+l) 

«=  ru(i-l,j)  -f  2(l-r)  u(i,j)  +  ru(i+l,j) 
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(3.3.4.10)  -rp2u(i-l,j+l)  +  2(l+rj32)  u(i,j+l)  -  r|32u(i+l,  j+l) 

=  r(32u(i-l,j)  +  2(l-rp2)  u(i,j)  +  r02u(i+l,j) 
where  u(i,j+l)  and  u(i,j)  are  not  adjacent  to  the  moving  boundary, 

-z  2 

The  local  error  in  both  of  the  above  equations  is  0((At)^)  +  0(At(Ax)  )  „ 
To  approximate  (3*3. 4. 6a),  (3.3.4.6b)  and  (3»3*4.7)j  u  is  expanded  in 
a  Taylor  series  about  the  point  (0,j+l)  in  both  x  and  t  directions 
and  equations  (3.3.4. l),  (3*3. 4. 2),  (3.3.4,6a),  (3,3.4.6b)  and  (3. 3. 4. 7) 
are  used.  We  obtain 

(3.3.4.11a)  (l+2r)  u(0, j+l )  -  2ru(l, j+l)  =  u(0, j)  +  2rAxp1f ( j+l)  , 

X(t)  =  0  , 

(3.3.4.11b)  (l+2rp2)  u(0,  j+l)  -  2rJ32u(l,j+l) 

=  u(0,j)  +  2rp2Axp2f(j+l),  X(t)  >  0  , 

(3.3.4.11c)  u(0, j+l)  =  g(j+l) 

where  the  error  term  in  (3.3.4,11)  is  the  same  as  that  for  (3«3«4.9)  and 
(3.3.4.10)  .  To  approximate  the  right  boundary  condition  we  suppose  grid 
points  to  be  at  (M+l,j)  and  (M+l,j+l)  where  M  =  B/Ax  and  suppose 
u  =  u(M“l,k)  for  k  =  j,  j+l  .  The  approximation  to  equation 

(3. 3* 4. 8)  is 

(3.3.4.12)  -ru(M-l,j+l)  +  (l+r )  u (M, j+l)  =  (l-r)  u(M,j)  +  ru(M-l,j) 


with  error  as  before. 


•  •  b  ■■■■  3  '  •  . 
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The  above  difference  equations  hold  only  for  those  points  not 
adjacent  to  the  moving  boundary.  Since  Ax  and  At  are  held  constant, 
the  moving  boundary  will  usually  lie  between  grid  points.  This  introduces 
complications  in  satisfying  equations  (3.3**+*9)  -  (3*3*^*12)  near  the 
boundary  curve. 

3.3*5  Difference  techniques  near  the  boundary. 

We  examine  the  difference  technique  which  Ehrlich  gives  for 
grid  points  adjacent  to  the  moving  boundary.  Assuming  we  have  an  estimate 
of  the  boundary  position  for  the  (j+l)st  row,  we  consider  four  cases 
each  of  which  are  handled  separately.  They  are; 

A.  boundary  curve  does  not  cross  a  space  grid  line, 

B.  boundary  curve  crosses  one  space  grid  line, 

C.  boundary  curve  crosses  two  space  gride  lines, 

D.  boundary  curve  crosses  three  or  more  space  grid  lines. 

i 

To  fix  the  notation  let  i  =  w  for  the  grid  point  on  the  moving  boundary 
or  to  the  left  of  the  boundary  on  the  (j+l)st  row.  If  the  boundary 
falls  on  a  grid  point,  the  corresponding  grid  line  is  not  crossed.  The 
difference  equations  take  care  of  the  case. 


Figure  3.3.5.! 


.?  >r*  "■  ?i'i  .  '  "  ■ 

■  "  *, .  \  . 


' 

' 
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Figure  (3.3*5«l)  is  an  example  of  case  A  where  s,  .  is  the  fractional 

j+i 

part  of  the  space  mesh  between  the  point  (w,j+l)  and  the  boundary.  The 
points  marked  with  *  require  special  equations.  Expanding  u(i,j)  in 
a  Taylor  series  about  (w,j+l)  in  both  the  space  direction  and  in  the  time 
direction  and  using  (3. 3*4. 2)  and  (3»3.4.3)>  we  obtain 


'2si+lrf3  2 

(3.3.5. 1)  1^g~r  - —  u(w-l,j+l)  +  (2rp  +sJ+1)  u(w,j+l) 


sj+iu(w’j)  +  uF 


j+l 


Similarly  on  the  right  of  the  boundary  we  obtain 

(30.5.2)  (2-sj+1)(l-2r-sj+1)  u(w+l,j+l)  -  2r(l-Sj+1)  u(w+2,j+l) 


(2-s  . 


j+l 


The  error  in  each  case  is 


)(!-Sj+i)  u(w+l,j)  +  2r 

0((At)^)-j-  O(AtAx)  . 


If  w  =  0,  the  boundary  lies  to  the  left  of  the  first  grid 
point.  Using  (3.3. 4. 2),  (3. 3.4.3)  and  (3.3.4.6b)  we  obtain 

(3.3.50)  (2r|32+s2j+1)  u(0,j+l)  =  s2+1u(0?j)  +  2r(32  +  2rp2s j+1^xp2f ( j+l ) 

which  replaces  equations  (3»3»5«l)  and  (3*3*5. 2)  . 


Figure  3. 3. 5. 2 


(  '■  >  ■  ,i v' 
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The  figure  above  illustrates  case  B  ,  where  s^  and  Sj+^ 
are  defined  as  before  and  p^  is  the  fractional  part  of  the  time  step 
between  (w,j+l)  and  the  moving  boundary „ 


(3. 3. 5. 4) 


Using  the  same  idea  as  in  case  A  we  obtain 


-rp2u(w-2, j+l)  +  (l+2r(32)  u(w-l,j+l)  -rp2u(w, j+l)  =  u(w-l,j) 


(3. 3. 5- 5) 


•2rg2  u(w-lt  j+l)  +  (2r@  +fj±I)  u<w>j+1> 


1+s 


j+l 


w 


(3.3.5 .6)  (2-sj+1)(l+2r-sj+1)  u(w+l,  j+l)  -  2r(l-sj+1)  u(w+2, j+l) 

=  (2“sj+i)(1_sj+i)  U(w+M)  +  2r  • 


Now  (3.3*5.6)  is  the  same  as  (3.3»5»2)  .  Equation  (3*3»5A)  is  a  back- 

2  2 

ward  difference  equation  with  error  0((At)  )  +  0(At(Ax)  )  .  Equations 

2 

(3.3.5. 5)  and  (3. 3. 5*6)  have  error  terms  0((At)  )  +  O(AtAx)  .  Assuming 
linearity  of  the  moving  boundary  between  two  successive  time  steps,  p^ 
can  be  written  as 


pw  = 


J±L 


4  0  . 


J±i  


=  1 “S  .  +s 


w 


j  j+l 


From  this  we  have 


. 


+ 

'  ■'  sS<  Q  +  { 
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If  Sj  =  I,  then  (3. 3. 5. 6)  becomes  (3.3.5. l)  with  u(w,j)  «  1.  We  note 
that  w  cannot  be  zero  for  this  case. 


Figure  3. 3. 5. 3 

The  difference  equations  for  case  C  ares 


(3. 3® 5®7) 

-r32u(w-3, j+1) 

(3.3. 5*8) 

”pw„lr32u(w-2,. 

(3»3«5«9) 

2 

-2r|3  s  .  . 

- L-J+L.  u(w 

1+sj+i 

(3.3.5.10) 

(2-sJ+1)(l+2r-i 

+  1 


I-i-s 


j+1 


’j+1' 


’  j+1' 


Equation  (3.3.5.IO)  Is  the  same  as  (3.3. 5. 2)  and  (3«3«5«6)  while  (3«3°5«1Q) 

is  the  same  as  (3«3«5«5)  with  a  modification  for  the  definition  of  p  as 

w 


In  (3. 3.5*9)  we  obtain, 


=  2-s.h-s 

P„  J  J+1 


w 


Figure  3. 3. 5. 4 

Case  D  is  treated  by  replacing  At  by  At/2  and  trying  to 
apply  one  of  cases  A,  B  or  C  . 

The  approximation  to  equation  (3»3»4.4)  is  now  considered.  This 
equation  locates  the  boundary  at  each  time  step  by  an  iterative  procedure. 
However, several  difference  schemes  are  necessary  because  of  difficulties 
encountered.  We  will  refer  to  the  figure  below  in  our  considerations. 


Figure  3-3. 5*5 
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To  approximate  the  space  derivatives  in  equation  (3. 3. 4.4),  a 
three  point  equation  is  used  on  each  side  of  the  moving  boundary,  the  three 
points  being  the  point,  and  the  next  two  points  (darkened  in  the  figure) „ 

We  obtain  the  following  approximations 


A, 


(3.3.5-11)  u  (x(t),t)  =  ^  +  0((Ax)2) 


where 


Ai  ■  u(w+i,j+i)  -  (A-J±l)u(w+2,j+i) 


l+2(l-s1+1) 


and 


(3.3.5.12) 


A*  ? 

ux(x(t)-0,t)  =  ~  +  Of  (Ax;  ) 


where 


s  .  ..  1+s  .  1 

A2  =  TIs±^“  u(w-l,j+l)  -  — J±-  u(w,J+l) 

j+1  j+1 


1+2s 1-t-l 

sj+l(1+“j+l) 


However,  certain  coefficients  in  (3«3«5.H)  and  (3. 3. 5. 12)  become  very 


large  as  approaches  0  or  1  „  This  necessitates  special  approxi¬ 

mations  as  the  moving  boundary  approaches  a  grid  point. 


Suppose  we  use  three  grid  points  on  each  side  of  the  bo.undary 
(squares  in  above  figure)  and  extrapolate.  The  difference  approximations 


are 


Yi.-i  >t  ■ 

■ '  .  .■ 
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a1 

(3.3.5.13)  ux(x(t)+0,t)  =  ^  +  0((Ax)2) 


where 


|[’(5“2sj+1)  u(w+l,j+l)  +  4(2-sj+1)  u(w+2, j+l) 


(3‘2sj+1)  u(w+3> j+l)l 


and 


(3.J.5.H) 


ux(X(t)-0,t)  =  ^  +  0((Ax)2) 


where 


A2  =  2^1+2sj+l^  u(w“2?j+1)  “  ^1+Sj+1)  u(w"l^j+l) 
+  (3+2s j+1)  u(w, j+l)]  . 


We  also  need  an  approximation  to  (3»3»^*J+)  at  the  start  of  melt¬ 
ing  (see  figure  below)  „ 


Figure  3.3. 5.6 


For  X(t)  >  0  ,  we  use  (3.3«^«6b)  „  If  we  expand  u  in  a  Taylor  series 

X 

about  (Q,t)  in  the  space  direction,  we  obtain 

A" 

ux(x(t)-0,t)  =  ~  +  O(AxAt)  +  0((Ax)2) 


(3.3.5.15) 
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where 


A£  =  -Axp  f(j+l)  +  (u(0,j+l)  -  u(0,j)) 

rp 


This  approximation  is  used  until  there  are  at  least  two  grid  points  (includ¬ 
ing  (0,t)  )  to  the  left  of  the  moving  boundary,  at  which  time  either  A^ 
or  A^  is  used.  It  is  necessary  to  use  a  special  equation  for  the  initial 
guess  on  the  time  step  immediately  following  the  start  of  melting.  This 
can  be  done  using  (3.3«^«6b)  and  A|  with  s  ^  =  1 

We  approximate  the  time  derivative  in  (3»3.^.*0  as  follows' 


(3.3.5.16) 


+  o(At)  . 


Since  X  is  the  distance  melted, 


X(j)  =  (Wj+Sj)  ^ 


where  VL  is  that  grid  point  on  the  moving  boundary  or  just  to  the  left 
of  it  at  the  jth  time  step.  Thus, 


(3.3.5.17) 


X(j+1)  -  X(j)  =  (Wj+1+S.+1)  Ax  -  (Wj+Sj)  Ax 


=  Ax  ((W.  ,-W.)  +  (s.  ,“S,))  =  Ax  (AW.  +As  .  ..  ) 

j+i  j  '  j+i  y  j+i  j+i' 


Thus,  from  (3.3.5.II),  (3. 3. 5. 12),  (3-3.5.16)  and  (3.3.5.17).  equat 
(3.3.^.^-)  is  approximated  by 


ion 


(3.3.5.18) 


=  AW  +As 
Ax  j+1  j+1 


At 


(Ax)' 


(A^-yAg)  +  0(.AxAt) 


‘ 


'  - .  * 
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Similar  equations  result  from  (3.3.5.13)  -  (3*3*5.15)  where  appropriate. 
We  note  that  0^  sj+l  ^  ^  an<*  ^j+1  an  ^nte8er* 


3.3.6  Method  of  solution 

Thus  (3. 3.^.1)  through  (3. 3. 4. 8)  are  now  replaced  by  (3»3.4.9) 
(3.3*j+»12)  and  (3*3.5.l)  -  (3.3*5»10)  where  appropriate  and  are  solved 
as  a  system  of  algebraic  equations.  Considering  (3. 3*^. 9)  *  (3.3»^«12) 
and  (3*3«5<>l)  -  (3 • 3 ® 5 • 10 )  one  notices  that  having  chosen  the  proper 
equations,  the  resulting  matrix  is  tridiagonal.  We  have  a  system  of  the 
form  for  each  j  . 


(5.3.6. 1) 


Bou(j,0)  +  CQu(j,l)  =  d(j,0) 

A.u(j,i-1)  +  BlU(j,l)  +  ClU(j,l+l)  =  d(j,i),  1  <  1  <  M-l 

A  u( j,M-l)  +  B  u(j,M)  =  d(j,M)  . 


If  one  carries  out  the  Gaussian  elimination,  eliminating  elements  below 
the  main  diagonal  of  the  matrix,  and  then  uses  back  substitution,  one  obtains 

u(j,M)  =  q(j,M) 


(3- 3.6.2) 


u(j,i)  =  q(j,i)  -  b(j,i)  u(j,i+l),  0  <  i  <  M-l 


where 


(3. 3. 6. 3) 


q(j,o) 


iLLOl 

B 


q(j,i)  = 


d(j,i)  -  Aiq(j,i-1) 


B 


±  -  A_Lb ( j , i“l)  ' 


1  <  1  <  M 
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and 


b(j,0)  = 


c 

o 


t 


o 


which  is  the  method  of  solving  for  the  (j+l)st  time  step  from  the  jth 
time  step. 

Ax,  At  and  M  are  chosen  and  then  (3*3*^»9)>  (3«3*^»Ha), 
(3.3.4.12)  and  the  initial  values  are  used  to  compute  u(j,i)  for 
0  <  j  <  M  and  i  >  1  using  the  above  technique  until  for  some  i  =  J, 
u(0,J+l)  >  1  and  u(0,J)  <  1  .  At  that  time 


r  or 


is  calculated  and  a  new  (j+l)st  time  step  is  computed.  This  becomes 
the  base  line  for  phase  II  . 

Setting  r"  =  r~r',  the  (j+2)nd  time  step  is  computed.  As 
a  first  approximation  to  the  slope  we  compute 


[|-(-3+4u(1,J+2)  -  u(2, J+2) )  +  yAxp2f  (j+2)  J  . 


Using  this  as  the  value  of  X(j+2)  =  W  +s  0  ,  the  (j+2)nd  line 

J-H& 

is  recomputed  using  the  appropriate  equations.  For  the  next  approximation 
to  X  the  following  criterion  is  used  where  j+1  replaces  J+2  for  the 
general  cases 


if  (Wj+1+Sj+1)  <  l  “se  in  (3.3.5. 18)  , 


:  .Ml 

’  ■■  .  '  .  % 
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if  \  <  «j+1  <  1  ,  use  A«  in  (3.3.5.18)  , 

if  0  <  sj+1  <  |  ,  Wj+1  >  1  ,  use  A^  in  (3.3.5*18)  , 

otherwise  use  (3»3»5.l8)  as  is  . 

For  the  (j+2)nd  time  step,  r"  replaces  r  in  (3*3*5- 18) » 


If  we  define 


k 

Zj+1 


>1 


t 


k  being  the  maximum  number  of  iterations,  (3.3* 6. 5)  gives  X(j+2),  then 
the  next  approximation  to  X  is 

Xk(j+1)  ■-  \  (Xk_1(j+1)  + 

for  the  general  case.  With  this  new  approximation,  the  line  is  recomputed 
using  the  appropriate  equations.  This  process  is  repeated  either  a  fixed 
number  of  times  or  until  convergence  to  some  criterion  is  reached.  Then 
the  (j+2)nd  step  is  assumed  known. 


Replacing  r"  by  r  ,  the  (j+3)rd  row  is  computed  by  first 
projecting  X(t)  linearly  from  X(j+2)  to  the  (j-j-3)rd  row  using 
a  slope  of  —4^)  as  computed  from  (3. 3.5«l8)  calling  this  X°(j+3), 
and  then  repeating  the  above  procedure.  Computation  continues  until 
t  >  t  ^  ,  a  previously  fixed  parameter. 

The  above  method  can  be  applied  to  problems  with  melting  tern*, 
perature  other  than  1  ,  The  method  appears  to  be  one  which  can  be 


■  i"  .  V  . 
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programmed  since  the  criteria  used  are  tested  readily  on  a  computer.  The 
problem  lies  in  choosing  appropriate  intervals  for  the  time  and  space 
variables . 

Ehrlich  gives  the  results  for  a  particular  problem  for  several 

values  of  Ax  and  At  for  which  — is  constant.  He  points  out 

(Ax)2 

that  one  is  usually  more  interested  in  the  accuracy  of  the  location  of 
the  moving  boundary  than  in  the  overall  solution.  He  finds  that  the 
method  give  smooth  convergence  to  the  solution.  For  a  fixed  mesh  size 
and  a  given  tolerance,  the  scheme  requires  a  varying  number  of  iterations 
depending  on  what  time  step  is  being  computed,  the  first  one  usually 
requiring  the  most  iterations,  and  which  one  of  (3»3«5«H)  -  (3®3«5»15) 
is  being  used.  He  finds  that  if  a  fixed  number  of  iterations  are  performed, 
then  the  locations  of  the  moving  boundary  for  varying  mesh  size  cannot  be 
distinguished,  Ehrlich  finds  his  method  stable  against  round-off  error 
since  the  round-off  error  grows  linearly  but  at  a  very  slow  rate  so  as 
not  to  affect  the  solution.  His  results  are  consistent  with  those  obtained 
by  Douglas  and  Gallie  (1955) •  They  do  some  numerical  integration  and  it 
is  difficult  to  compare  the  results  since  it  is  not  known  how  much  of  the 
error  is  due  to  numerical  integration, 

3.3.7  Summary 

The  main  difficulty  in  finding  numerical  solutions  for  the  moving 
boundary  problem  is  finding  a  suitable  algorithm  for  locating  the  boundary. 
The  non-analytic  techniques  are  A.  E.  approximations  to  the  D,  E, 
and  the  initial  and  boundary  conditions.  It  appears  that  Ehrlich's 


i.j'  x.  •.■'.v,:: :  .  ->3*  >Vs.2iJri  X  Xx;.wx  -;x  .xx1:’ 
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method  is  most  useful  since  it  gives  an  algorithm  which  is  machine 
computable  and  gives  results  comparable  to  those  obtained  by  other  more 
complicated  techniques. 
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CHAPTER  4 

CONCLUSION 


In  the  preceding  chapters,  several  proofs  of  existence  and 
unicity  of  the  solution  of  a  moving  boundary  problem  have  been  outlined. 
The  presentation  of  a  formal  solution  again  demonstrates  the  existence 
and  unicity  of  a  solution  to  the  problem.  The  formal  solution  yields  a 
functional  equation  for  the  boundary. 

Several  analytic  and  non-analytic  approximations  to  solutions 
have  been  examined.  In  the  process  of  finding  analytic  approximations, 
the  repetitive  evaluation  of  coefficients  of  a  power  series  for  which 
no  radius  of  convergence  has  been  given  proves  to  be  tedious.  The  main 
difficulty  in  both  analytic  and  non-analytic  procedures  is  locating  the 
moving  boundary.  Non-analytic  techniques  used  in  linear  problems  may  be 
employed  except  near  the  moving  boundary  where  non-linearity  of  the 
problem  requires  special  techniques.  The  usual  procedure  is  to  estimate 
the  location  of  the  boundary  by  trial  means  and  then  to  use  an  iterative 
scheme  to  converge  on  the  boundary. 


V  ■  '  i  r.i.hu 
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